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This paper describes a method for determining the dynamical interaction between 
extended halo and spheroid components and an environmental disturbance. One finds 
that resonant interaction between a galaxy and passing interlopers or satellite com- 
panions can carry the disturbance inward, deep inside the halo, where it can perturb 
the disk. 

Applied to the Milky Way for example, the LMC and SMC appear to be sufficient 
to cause the observed Galactic warp and possibly seed other asymmetries. This is a 
multi-scale interaction in which the halo wake has a feature at roughly half the satellite 
orbital radius due to a 2:1 orbital resonance. The rotating disturbance then excites 
an m = 1 vertical disk mode which has the classic integral-sign morphology. A polar 
satellite orbit produces the largest warp and therefore the inferred LMC orbit is nearly 
optimal for maximum warp production. 

Both the magnitude and morphology of the response depend on the details of 
the disk and halo models. Most critically, a change in the halo profile will shift the 
resonant frequencies and response location and consequently alter the coupling to the 
bending disk. Increasing the halo support relative to the disk, a sub-maximal disk 
model, decreases the warp amplitude. 

Finally, the results and prognosis for N-body simulations are discussed. Discrete- 
ness noise in the halo, similar to that due to a population of lO^M© black holes, can 
produce observable warping. Details are discussed in an associated paper. 

Key words: Galaxy: halo, structure — galaxies: halos, kinematics and dynamics — 
Magellanic Clouds 



1 INTRODUCTION 

Even casual examination shows that most disk galaxies are 
not truly symmetric but exhibit a variety of morphological 
peculiarities of which spiral arms and bars are the most pro- 
nounced. After decades of effort, we know that these features 
may be driven by environmental disturbance acting directly 
on the disk, in addition to self-excitation of a local distur- 
bance (e.g. by swing amplification, Toomre 1981, Sellwood 
& Carlberg 1984). However, all disks are embedded within 
halos and therefore are not dynamically independent and 
will respond to asymmetries and distortions in the halo, as 
well. 

Until recently, conventional wisdom was that halos 
acted to stabilize disks but otherwise remained relatively in- 
ert. The argument behind this assumption is as follows. Ha- 
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los, spheroids and bulges are supported against their own 
gravity by the random motion of their stars — a so-called 
"hot" distribution (e.g. Binney & Tremaine 1987). On all but 
the largest scales, they look like nearly homogeneous thermal 
baths of stars. Because all self-sustaining patterns or waves 
in a homogeneous universe of stars with a Maxwellian veloc- 
ity distribution are predicted to damp quickly (e.g. Ikeuchi, 
Nakamura & Takahara 1974), one expects that any pattern 
will be strongly damped in halos and spheroids as well. How- 
ever, recent work suggests that halos do respond to tidal 
encounters by companions or cluster members and are sus- 
ceptible to induction of long-lived modes due to their inho- 
mogeneity. These modes are at the largest scales for which 
self gravity is most effective. In particular, if halos are large 
as many currently estimate, halo-halo interactions in groups 
or clusters will be frequent and much more common than 
disk-disk interactions. Because non-local coupling can carry 
a disturbance to small radii, transient halo structure may 
be sufficient to trigger disk structure even if the halo pat- 
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tern subsequently damps away. The non-locality of the re- 
sponse was predicted by Weinberg (1986, 1989) and verified 
in n-body simulations by Hernquist & Weinberg (1989), and 
Prugniel & Combes (1992). Application to self-gravitating 
fluctuations is described in an associated paper (Weinberg 
1997). 

Similarly, we expect that a companion can continuously 
re-excite structure. The orbit of the companion decays by 
dynamical friction. The response of the dark halo to the in- 
terloping satellite can be thought of as a gravitational wake 
(e.g. Mulder 1983); since the wake trails and has mass, it ex- 
erts a backward pull on the satellite. This view of dynamical 
friction reproduces the standard approach but shows that 
the halo responds with structure whose mass comparable 
to the satellite. This general approach has been tested and 
compared with n-body simulations in a variety of contexts 
with good agreement. Recent work (Weinberg 1995b, Paper 
I) suggests that the Magellanic Clouds use this mechanism 
to produce distortions in the Galactic disk sufficient to ac- 
count for both the radial location, position angle and sign of 
the HI warp and observed anomalies in stellar kinematics. In 
clusters, this mechanism is mostly likely the culprit behind 
galaxy harassment (Moore et al. 1996). 

Here, we develop this suggestion and present a formal- 
ism for exploring the dynamics of low-amplitude interac- 
tions that can lead to significant long-term evolution. As 
an example throughout, I will focus on disk warping and 
on the Milky Way — LMC interaction presented in Paper I. 
Although there are some general trends, a relatively large 
amplitude disk response tends to be the result of a conspir- 
acy between frequencies. Rather than present an exhaustive 
parameter survey, we explore some simple scenarios. Even 
though the warp amplitudes vary with the details of the 
galaxy profiles, astronomically interesting amplitudes are 
generally produced. Sections ^ and |^ describe the galaxian 
models and the method. We will use a numerical perturba- 
tion theory which is well suited to describing weak coherent 
perturbations. The main results are in ^ which begins by 
examining an example distortion of the halo by a compan- 
ion satellite and traces its influence on both the disk warp 
and in-plane disk distortions. Readers may skip the technical 
detail without loss of continuity by turning to §^ after the 
introduction to §sec:method. The results section is followed 
by a discussion of the range of effects using other models 
and rough generalizations (§^. It is tempting and desirable 
to follow this up with n-body simulation but simulations of 
weak distortions with large dynamic range is a significant 
challenge and will require very large particle numbers to re- 
cover signal (§h). We end with a summary and outline for 
future work (§m). 



2 GALAXY MODELS 

In order to explore interactions between components, we 
need self-consistent multi-component galaxy models. Fully 
self-consistent disk-in-halo models are generally made pre- 
scriptively rather than constructively because of the diffi- 
culty in determining distribution functions. For example, 
a regular disk profile embedded in an a regular spherical 
halo will generally not yield a self-consistent regular sys- 
tem demanding computationally intensive techniques such 



as Schwarzschild's method (1979). Most often in the litera- 
ture, the Jeans' moment equations are used to construct an 
n-body disk in a given halo or spheroid and the resulting 
distribution is allowed to phase mix to equilibrium. In this 
paper, we adopt a hybrid approach suited to both the an- 
alytic perturbation theory described in ^ and the n-body 
simulations described in ^ The prescription is as follows: 

(i) Choose a spherical halo and two- or three-dimensional 
disk profile. 

(ii) Assume that the disk does not affect the halo pro- 
file and construct a disk phase-space distribution function 
for the disk in the halo. This would rigorously obtain for 
Mhaio 3> Mdiak- The disk distribution function is com- 
puted by a quadratic programming scheme similar to that 
discussed by Dejonghe (1989, see Appendix ^). The term 
halo here includes the entire hot component: bulge, stellar 
spheroid and dark matter. The halo distribution function is 
assumed to be known. If it is not, it may be constructed 
using any of the established techniques (e.g. Eddington in- 
version, generalized integral inversion, atlas method, etc.), 
ignoring the disk. 

(iii) The approximate distribution functions for each com- 
ponent are now known. For an n-body simulation, these may 
be directly realized by a Monte Carlo procedure. 

The realized phase space distribution will not be in strict 
equilibrium. However, as long as the force of the halo domi- 
nates the disk at large radii, and the disk dominates its own 
gravity inside its scale length for realistic parameters, the 
initial construction is mostly likely close to equilibrium. Nu- 
merical experiments support this conjecture. Nevertheless, a 
mild inconsistency is of minor consequence for t he numerical 
method presented in §^ (see discussion in §3.2). 

For analytic convenience, we use King models of varying 
concentration and scale to represent the dark halo, a Hern- 
quist model to represent the bulge/spheroid, and Hunter 
(1963) polynomial disk profiles. The latter choice is moti- 
vated by the numerical analysis required to compute the ver- 
tical disk response following Hunter & Toomre (1969, here- 
after HT). Hunter's polynomial disk models are less centrally 
concentrated and fall of more steeply then the standard ex- 
ponential disk. Hunter & Toomre modified these disks to 
better fit observed profiles by subtracting off a low-order 
contribution; they denoted these models with the suffix "X" . 
The Hunter- Toomre 16X model is fair approximation to the 
exponential disk with a scale length of Rmax/G; this cor- 
responds to a = 3.5,4.5kpc for Rmax = 21,28kpc. These 
profiles are compared in Figure |l|. For the n-body tests de- 
scribed in §^ we used a rigid bulge component to stabilize 
the inner disk (which is otherwise bar unstable). 

The extent of Milky Way disk remains under debate. 
Robin, Creze & Mohan (1993) present evidence for an edge 
to the stellar disk at 14 kpc. This is consistent with esti- 
mates based on molecular tracers (Wouterloot et al. 1990, 
Digel 1991, Heyer et al. 1997) which imply an edge of about 
14 kpc or somewhat greater. The outer atomic gas matches 
on to the inner component continues to 30 kpc at signif- 
icant surface densities (Kulkarni, Heiles & Blitz 1982). A 
more recent analysis (Diplas & Savage 1991 ) reports an 
exponentially distributed neutral hydrogen disk out to at 
least 30 kpc. We will adopt the Rmax ~ 28 kpc model be- 
cause it better represents the global extent of the disk even 
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Figure 1. Left-hand (right-hand) panel compares the exponential disk to the Hunter-Toomre 16X model with identical mass inside of 
-R = 3 (-R = 4) corresponding to approximately 21kpc (28kpc). 



though the scale length required is a bit large. In calcula- 
tions below, the disk model is assigned M — 1 and this 
corresponds to the estimated 6.0 x 10^° M0 disk (Binney 
& Tremaine 1987). The radial scaling is one unit for each 
7kpc. The satellite orbit is defined by its energy and an- 
gular momentum in the halo model. For each halo model, 
the orbit is assigned a pericenter at 50 kpc and apocenter at 
100 kpc. The mass of the LMC in these units is 0.25 (0.1) 
for Mlmc = 1.5 X 10^° Mq (6 x 10^ Mq) (Meatheringham 
1988, Schommer et al. 1992). For this study, we lump the 
Small Cloud together with the LMC and do not consider 
the possibility of a distinct SMC orbit. In all that follows, 
we will scale results to the Milky Way. 

Given the disk profile, the King radius and mass are 
constrained by the demand for a flat rotation curve, at least 
for R < 50 kpc. The choice of a Wo ~ 3 model with trun- 
cation radius at 200 kpc and mass of 10 disk masses results 
in a Milky Way mass of 3.3 x 10" inside of 50 kpc. This 
is a bit smaller than but comparable to current values; e.g. 
Kochanek (1996) estimates 4.9 x lO" Mq. We also wiU con- 
sider halo mass ratios of 15 and 20 disk masses. The rotation 
curves remain approximately flat for both cases and give 
masses within 50 kpc of 4.7 x lO" M© and 6.0 x lO" Mq, 
respectively. Finally, the timing argument and its recent gen- 
eralizations (Peebles 1995 and references therein) estimate 
a Milky Way mass of 2 x 10^^ M0 . This is a constraint not 
a target for our models since one can add mass beyond 100 
kpc without affecting any of the dynamical arguments con- 
sidered here. Our most massive model is below this limit 
with a mass of 1.3 X 10^2 Mq. 

The perturbing satellite orbit is chosen to match the 
LMC orbit inferred by Lin et al. (1995) and assigned an 
orbit consistent with the chosen halo profile and scaling 
given in the previous section. The clouds are assumed to be 
near pericenter with Rlmc ~ 49.5 kpc. The space velocity 



inferred from proper motion, radial velocity measurement, 
and a distance estimate immediately determines the orbital 
plane as described in Paper I. In a Galactocentric coordi- 
nate system with the LSR along the a;-axis {x = —8.5 kpc) 
and moving towards positive y, one may straightforwardly 
derive the following instantaneous orbit: —76 ± 13° inclina- 
tion, —82 ±10° longitude of ascending node, —36 ±3° argu- 
ment of perigalacticon. However, considered in these Galac- 
tic coordinates, the Milky Way disk rotates clockwise and 
has negative angular momentum. In the development below, 
disks are assigned positive angular momentum with counter- 
clockwise rotation. Our models may be transformed to the 
Milky Way system by reflection through the y-z plane. 

We will adopt a standard galaxy model with a King 
model halo with tidal radius Rt = 28 (200 kpc), log c = 0.67, 
and a mass of 10 times the disk mass; this is "maximal" disk 
model. The bulge is Hernquist model with scale length 0.2 
(1.4 kpc) and mass of 20% of the disk (1.2 x 10^° Mq). The 
disk is the Hunter-Toomre 16X model with Rmax = 4 (28 
kpc). In this system, one velocity unit is 350kms~^. The 
rotation curve for this model rises from 0.6 at i? = 0.5 to 
0.7 mt R = 1.8, drops slowly to 0.6 at = 10 (70 kpc) and 
drops off more rapidly beyond this point. Although better 
fits to the observed Milky Way rotation curve are available, 
our goal of understanding the underlying mechanism and 
the computational simplicity of these components supports 
our choice. Finally, the standard model includes a satellite 
with an LMC orbit. The magnitude of the response scales 
with satellite mass and need not be chosen a priori. 



3 A FORMALISM FOR MULTI-SCALE 
INTERACTIONS 

A full treatment requires the dynamical coupling of the mul- 
tiple time scales and multiple length scales of the external 
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disturbance and the galaxian components discussed above. 
Relevant characteristic length and time scales may differ by 
an order of magnitude between satellite and halo or disk or- 
bits. In addition, we will see that the halo disturbance may 
be relatively weak and a significant perturbation of the outer 
disk at the same time. These multiple-scale weak regimes are 
a challenging task for an n-body computation. However, this 
class of problems is ideally suited to linear techniques and 
the work here will use the expansion technique known as 
the matrix method. Although the matrix method is compu- 
tationally intensive, it is no more so than n-body methods 
and is practical on current workstations. In this section, I 
will give a brief overview of the general method with details 
on posing and implementing the coupled response solutions 
in the references cited below and in the Appendix. 

In short, the matrix method represents the response 
of a galaxy to an external perturbation by a truncated se- 
ries of orthogonal functions, similar to those one would use 
to solve an electrostatics problem. The perturbation is also 
represented by this series and the temporal dependence of 
each coefficient is Fourier transformed to a (complex) fre- 
quency distribution. The response of the galaxy to one of 
the orthogonal functions at a particular forcing frequency is 
then computed in the continuum limit using the coUisionless 
Boltzmann (Vlasov) equation. The entire procedure is anal- 
ogous to signal processing in Fourier space. Pursuing the 
analogy, we now do the inverse transform. The response to 
any perturbation, the weighted superposition of the response 
to each basis function, is then a matrix equation. Finally, to 
get the full time dependence of the response, one resums the 
solutions to the matrix equation at each frequency weighted 
by the Fourier coefficients. 

This method assumes that the perturbation is small 
enough that the overall change to the structure of the galaxy 
is small. In this limit, the method has the advantage of ac- 
curacy and sensitivity to the large scale structures of in- 
terest. For contrast, the n-body simulation determines the 
response of a galaxy to a perturbation by solving the equa- 
tions of motion for a representative set of orbits. The orbit 
is advanced in a fixed potential for a short time interval and 
the gravitational potential or force is then recomputed. The 
simulation works well for large perturbations but because 
the simulation uses a finite number of particles, fluctuation 
noise limits the sensitivity to small amplitude distortions. 
The matrix method nicely complements the n-body simula- 
tions, excelling in the regimes where the n-body simulation 
are suspect. 

Historically, the approach is related to the treatment of 
general eigenvalue problems described in the mathematical 
physics literature (e.g. Courant & Hilbert 1953, Chap V). 
The matrix method in stellar dynamics had varied applica- 
tions beginning with Kalnajs (1977) who investigated the 
unstable modes of stellar disks. Polyachenko & Shukhman 
(1981) adapted the method to study a spherical system (see 
also Fridman & Polyachenko 1984, Appendix) and it was 
later employed by both Palmer & Papaloizou (1987) in the 
study of the radial orbit instability and by Berlin & Pego- 
raro (1989) to study the instability of a family of models 
proposed by Berlin & Stiavelli (1984). In addition to Paper 
I, Weinberg (1989, Paper II) used the matrix formulation to 
study the response of a spherical galaxy to an encounter with 



a dwarf companion and Saha (1991) and Weinberg (1991) 
investigated the stability of anisotropic galaxian models. 

3.1 Mathematical overview 

The response of a galaxy initially in equilibrium to a grav- 
itational interaction with a companion is described by the 
simultaneous solution of the Boltzmann and Poisson equa- 
tions. The simultaneous system is a set of coupled partial 
integro-differential equations. However if the orbits in each 
component are regular, any phase-space quantity — such as 
density and gravitational potential — may be expanded in 
a Fourier series in the orbital frequencies. Truncating this 
expansion, the quantity may be represented as a vector 
of Fourier coefBcients; this is standard practice in flltering 
and approximation theory and canonical perturbation the- 
ory (e.g. Lichtenberg & Lieberman 1983). In Fourier space, 
the Boltzmann PDE becomes an algebraic integral equa- 
tion. The system is further simplifled if the basis functions 
are chosen to satisfy the Poisson equation explicitly. After 
a Laplace transform in time, the remaining solution of the 
Boltzmann equation becomes the solution of a matrix equa- 
tion, each column describing the response to a particular 
basis function. See references cited abouve for mathemati- 
cal detail. 

To give the flavor of its use, consider two interacting 
galaxies. Denote the expansion of the perturbation potential 
caused by the companion galaxy as vector b. Then the direct 
response of the galaxy, vector a may be written: 



a = R b. 



(1) 



The matrix R, the response operator, implicitly contains the 
time-dependence of the perturbation as well as the dynamics 
of the Boltzmann equation. If we are interested in the self- 
gravitating response to the perturbation, we need to solve: 



a = R ■ (a + b) . 



(2) 



In words, equation (^j states that the self-consistent reaction 
of a galaxy to a perturbation is the gravitational response 
to both the perturbing force and the force of the response 
itself. Equations (^ and (^) result from the Laplace trans- 
form of the Boltzmann equation and therefore represent a 
particular frequency component, a = a(s). Therefore, the 
solution of equations (|l|) or (^ requires an inverse Laplace 
transform to recover its explicit time dependence (see Paper 
II for details). See Nelson & Tremaine (1997) for a general 
discussion of the response operator. 

3.2 Combining multiple components 

This approach is easily extended to flnd the simultaneous 
self-consistent response of several galaxian components. For 
example, let us consider a halo and a disk; any number of 
components may be combined similarly. Each component's 
distribution function solves a Boltzmann equation and is 
coupled to the others only through the total gravitational 
potential. If the interaction between the components is arti- 
flcially suppressed, the simultaneous solution is that of the 
following augmented matrix equation: 



Rh 








Rd 



+ 



(3) 
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where the subscripts h and d stand for the halo and disk, b 
is the external perturbation, and is defined to be the ma- 
trix with the same rank as R and all elements zero. Equa- 
tion (^) is simply two stacked versions of equation (^. The 
off-diagonal partitions in the augmented matrix, now 0, de- 
scribe the mutual interaction between components. 

To allow the components to interact, we project the 
halo response onto the disk basis and then add this to 
the right-hand-side of the disk response equation and vice 
versa for the response of the halo to the disk. Letting the 
matrices that perform these projections be Phd and Pdh, we 
may write the fully coupled set as 



( 











Rd 










J^d 


R 


h 


RdPdh 


Rji 








Rd 
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Phd 


Pdh 


1 


hh\ 




bd 





+ 



Rd 



( Kh I \ fhh\ 
V I Rd J '{hdj 



(4) 



The first term on the right-hand-side describes the mutual 
self-gravitating response due to each component and the sec- 
ond term describes the response of each component to the 
external perturbation. 

We may straightforwardly isolate effects of interest by 
coupling or uncoupling components or by including or sup- 
pressing self-gravity. For example, we may consider the re- 
sponse of the disk to a halo wake, but without the back 
reaction of the halo to the disk by setting the upper right 
term in the augmented response matrix to zero. Or, if we 
want to limit the disk response to forcing by the halo alone, 
the direct response of the disk to the perturbation may be 
left out by setting bd to zero, yielding 



Rfe 





RdPdh 


Rd 



Rfe 








Rd 



.(5) 



This is equivalent to first solving a.h = R^ • (a^, + b^ ) and 
then ad = Rd • (ad + Pdh8Lh)- This formalism couples the 
components by their perturbation from equilibrium, not the 
fully gravitational attraction. For this reason, any mild de- 
viation from perfect self consistency produced by the pre- 
scription described in ^ does not cause any problems and is 
unlikely to be significant. Most of the computational work is 
in producing the matrices R. Afterward, all of the variants 
may be studied with little additional effort. 



3.3 Satellite perturbation 



In order to apply the technique described in §3.2 to a per- 
turbation by an orbiting satellite, we need to expand its 
potential in the chosen basis to get t he pe rturbation vec- 
tor b. This expansion is described in §3.3.1 and applied in 



3.3.2 



Perturbation by an interloping galaxy may be treated 



similarly but is not described here. 



3.3.1 Fourier expansion 

Our complete set of basis functions are pairs of functions, 
{pi,di), which solve Poisson's equation, V^di — AnGpi, and 
are biorthogonal: 



I drp*{r)dj{r) = Sij. 



(6) 



The potential for the arbitrary point mass may be expanded 
directly in a biorthogonal harmonic series: 



Im, 



^(r) = > ;y,™(e,0)^6i'"(t)prw 



(7) 



where 



bTit) 



d'rYUeA)p'r'{r)S'l 
= YU0it),4>{t))p'r'ir{t)) 



(8) 



where r(t) describes the orbit of the satellite. Alternatively, 
since fe'™ is an implicit function of time through r, we may 
expand in an action-angle series which makes the time de- 
pendence explicit. This gives 



E 



im^f il20L 



Y,^(7r/2,0)W^, 



lli* ^i{liQ.i-\-l2^2)'l^ 

I I2 TTL 



(9) 



where the coefficient Wij^^ depends only on the energy and 
angular momentum of the satellite orbit. Derivation of equa- 
tion (P) is given in the Appendix. 

The perturbation vectors b in equation (^) are given by 
the Laplace transform of b'/^{t) from equation (Q) for each 
term in equation (^). The Laplace transform of bl'"{t) in 
this form is trivial. Because the physical measurable must be 
real, we can simplify the analytic computation for m 7^ by 
using only m > terms and adding the complex conjugate 
in the end. 

3.3.2 Application to response calculation 

Finally, recovery of the response, a, in the time domain re- 
quires the inverse Laplace transform of the following solution 
for a from equation fl): 









^) 
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Rd 
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RdPdh 


H-d )\ 
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(10) 



hh 
^bd 

or, in compact form, 
a''"(s) = ©~'(s)■7^(s)■b''" 

The matrix in large square brackets in this equation is the 
dispersion relation; its determinant vanishes at eigenmodes. 
The assumption that our multi-component galaxy is stable 
ensures that the inverse of this matrix has no poles in the 
complex half plane with K(s) > 0; poles on the half plane 
with K(s) < correspond to damped modes. Elements of the 
second term matrix have at least one pole on the imaginary 
s-axis due to the harmonic forcing by the satellite. In the 
end, the inverse transform may be simply evaluated by de- 
forming the integration contour through the imaginary axis 
and taking the time-asymptotic limit (see §^for details). Af- 
ter many satellite orbits, the harmonic-forcing contribution 
dominates all but a very weakly damped mode. 
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3.4 Disk bending 

A diflerential vertical force applied by the external pertur- 
bation and the halo wake can warp the disk plane. HT de- 
scribe a linearized solution for the dynamical evolution of 
the plane for an isolated disk. The bending analysis here 
uses the formalism developed by HT with several modifi- 
cations. First, because our equilibrium disk models are em- 
bedded in an external halo, we must retain their equation 
(12) rather than simplify using relationships based on the 
specific form of the background model. Second, we solve the 
linearized equations of motion under a forced disturbance 
(their eq. 19) by L aplace transform for consistency with the 
approach in §3.2, The vertical force follows directly from 
the expansion coefficients describing the external perturba- 
tion, equation (^ and the halo wake, equation (^^. The 
back-reaction to the in-plane distortion is included and is 
a relatively minor contribution to the total response. This 
calculation does not consider the back-reaction of the halo 
to the vertical distortion (Nelson & Tremaine 1995). 

For reasons described in HT, their polynomial disk mod- 
els are well suited to numerical analysis and adopted here as 
described in §p. I also tried exponential disks with different 
basis sets but could not find an alternative which allowed an 
accurate computation of the height alone, rather than the 
combining height times the surface density. 



4 WAKES AND WARPS 
4.1 Wake in halo 

There are two contributions to the disk disturbance: the di- 
rect tidal force of the satellite on the disk and the force of 
disturbance excited by the satellite in the halo — the halo 
wake — on the disk. The strength of the halo wake is pro- 
portional to the mass of the satellite and mass density of 
the halo. However it depends critically on the orbital struc- 
ture of the halo because a particle's orbit will respond most 
strongly near resonances between the satellite's and halo 
particle's orbital frequencies. Co-orbiting trajectories will 
have the strongest response but the total mass involved is 
small for the standard model. Higher-order resonances will 
be weaker but occur at smaller galactocentric radii where 
the mass density is high. The wake is the product of these 
competing efi'ects and, generally, the wake peaks far inside 
the satellite orbit. 

As an example. Figures ^ and |^ show the space density 
distortions induced by the LMC orbit in the standard King 
model halo in the orbital plane. The satellite has pericenter 
at f? = 7 and apocenter at i? = 14 and here, orbits in the 
counter-clockwise direction. Pericenter is X = 7, y = 0. If 
the satellite were completely outside the halo, the dipole re- 
sponse would be a linear displacement representing the new 
center of mass position. The wake would be proportional to 
—dp{r)/dr ■ e where e is the unit vector from the halo to 
the satellite center (Weinberg 1989). In Figure |2|, the satel- 
lite is inside the halo, and the wake deviates from the pure 
displacement. The amplitude of the quadrupole (Fig. ^) is a 
factor of roughly two smaller than the dipole. The dominant 
wake is near the satellite pericenter as expected, but note 
the inner lobe of the wake at roughly half the pericenter 
distance. This is due primarily to the 2:1 resonance between 
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Figure 2. Wake in the fiducial halo due to orbiting satellite for 
I = m = 1 (dipole). Nine linear spaced contours of overdensity 
(solid) and underdensity (dotted) shown. Each successive panel 
shows the wake radial phase as labeled, with pericenter at phase 
and apocenter at phase tt. Satellite locations at the four phases 
are {X,Y) = (7.0, 0.0), (0.2, 11.6), (-8.5, 11.2), (-11.3, 2.9). The 
outer and inner circles indicate the pericentric and half-pericentric 
radii, respectively. 



satellite and halo orbital azimuthal frequencies. Although 
the relative density of this inner lobe is smaller than the 
primary outer one at pericenter (phase 0), both the prox- 
imity and spatial structure causes the force from the inner 
wake to dominate over direct force from the satellite. As 
the satellite approaches pericenter (phase 37r/2), these in- 
ner lobes become relatively stronger and can dominate the 
response. The inner galaxy wake is weaker past pericenter 
(phase 7r/2). 



4.2 Vertical force on disk due to wake 

In the absence of the halo, the first multipole contribut- 
ing to the differential or tidal acceleration of the disk is at 
quadrupole (Z = 2) order. It is straightforward to convince 
oneself of this fact: the I — term is constant yielding no 
force, the 1 = 1 term is proportional to r yielding a spatially 
constant force, and therefore the 1 = 2 term provides the 
lowest order differential force. Because the warp has m = 1 
symmetry, the dominant warp inducing term will be Z = 2, 
jmj = 1. Similar symmetries apply for the action of the per- 
turbed halo on the disk. Including the halo warp, the dipole 
still can not produce the classic odd integral-sign warp but 
causes an even distortion. The lowest order halo wake that 
can excite an integral-sign warp is also the I = 2, \m\ = 1 
component. To illustrate the three-dimensional structure. 
Figure ^ renders the isosurface corresponding to 75% of peak 
amplitude wake contoured in Figure M. The wake is symmet- 
ric about the satellite's orbital plane and this plane is easily 
visualized. A wake must asymmetric about the z axis to 
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Figure 3. As in Fig. g but for I = m = 2 (quadrupole)halo 




Figure 4. A three dimensional rendering of the I = \m\ = 2 halo 
wake shown in Fig. ^). The isolevel shown at 75% of the peak 
amplitude. Overdensity and underdensity are shaded light grey 
and dark grey, respectively. The wire-frame outline extends ±10 
units in x, y, and z and the x-y-z axes are shown with the z- 
direction along the vertical. The satellite pericenter is a bit over 
7 units and shown as a small sphere. 

cause a differential vertical acceleration of the disk; in other 
words, a satellite in the disk plane produces no warp. The 
maximum vertical force occurs when the pattern shown in 
Figures and ^ is oriented perpendicular to the disk plane 
(see Fig.^ . The orientation of the LMC disk plane is nearly 
ideal for producing a disk warp. 



Figure 5. Warp amplitude as a function of orbital plane incli- 
nation. A polar orbit produces the maximum warp with heights 
473 pc (1181 pc) for LMC mass of 6 X lO^ Mq (1.5 X IQio Mq). 
The filled square shows the elevation of inferred LMC orbit. 




Figure 6. Vertical force on disk plane due to orbiting satellite 
only for harmonics I = 2, |m| = 1. Shown midway between peri- 
center and apocenter, at radial phase 7r/2 (cf. Fig. 

Figures describes the vertical force on the disk plane 
due to both the satellite alone and the halo response to the 
satellite for the standard model. Although the net force from 
the satellite is similar in magnitude to that from the halo, the 
satellite force is varies linearly with distance and only gives 
rise to a uniforming tilting of the disk. However, the spatial 
structure in the halo force causes a differential bending of 
the disk with vertical m = 1 symmetry: a warp. 

4.3 Vertical response of the disk 

One can think of the vertical response as a superposition or 
packet of modes. Bending modes in the absence of a halo 
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Figure 7. As in Fig. ^ but for the halo response to the satellite. 
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Figure 8. Tipping modes for the Hunter A'^ = 16 disk with edge 
at ij = 4 (approximately 28 kpc) embedded in a Wo = 3 King 
model with Rt = 28 (approximately 200 kpc scaled to the Milky 
Way). Modes are labeled in the notation of HT. 



a; 2 - 




Figure 9. Evolution of the tipping modes for the Hunter- Toomre 
16X disk embedded in a Wo = 3 King model with Rt = 28 
(approximately 200 kpc scaled to the Milky Way) as a function 
of halo mass. Curves are labeled by halo-to-disk mass ratio. At 
zero halo mass, the mode is bodily tipping of the disk at zero 
frequency. Each curve parallel to the y axis describes the mode 
with an increasing halo mass fraction shown on the x axis, up to 
a possible total of M = 30 times the disk mass. 
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have been described by HT. The addition of the halo shifts 
and slightly modifies the shape of these modes although they 
are qualitatively similar to those described in Figure 3 of 
HT. The corresponding modes in the standard halo model 
{Wo = 3 King model halo with truncation radii at Rt = 28 
[200 kpc]) are shown in Figure |^. In the absence of the halo, 
the zero frequency m = 1 mode is a bodily tipping of the 
disk. As the halo mass increase relative to the disk mass, the 
shape of the modes change. The tipping mode, for example, 
evolves into a distortion with warp-like structure. For the 
standard (Hunter- Toomre 16X) disk, the tipping mode is the 
only discrete mode and appears to be a prominent feature of 
the response discussed below. Its shape for increasing halo 
to disk mass ratio is shown in Figure |^. The tipping mode 
is distorted from linear to an integral-sign shape. 



Figure 10. Vertical response of the Hunter- Toomre 16X disk 
due to the combined vertical force from the standard satellite 
and halo. Scaled to the LMC and Milky Way following ^, the 
peak height is 420 pc (1050 pc) for the lower and higher LMC 
mass estimates, respectively. 



The combination of the force exerted by the satellite 
and the satellite-induced halo wake excites a vertical re- 
sponse in the disk. A strong vertical disk response obtains 
for near commensurable halo-wake pattern speeds and disk 
bending mode frequencies. These are accidental resonances 
in the sense that their existence is circumstantial and not 



the result of a tuning mechanism. Figure 
sponse to the combined force from Figures 



10 shows the re- 
tI The vertical 
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Figure 11. m = 1 density distortion in the disk plane due to the 
satellite and halo perturbation. Contours are evenly spaced from 
10% to 90% of maximum. 

distortion of the midplane within a galactocentric radius of 
roughly 10 kpc is less than 100 pc, however, the warp reaches 
a peak height of about 400 pc (1.0 kpc) at Rg ~ 20 kpc for 
LMC mass of 6 X 10*^ Mq (1.5 X 10^° Mq). In many cases, the 
warp has a local maximum in the outer disk, corresponding 
to the location maximum curvature in the integral sign. At 
larger radii the mode may turn over and reach a global max- 
imum amplitude at the very edge of the disk, corresponding 
to the ends points of the integral sign. We will refer to this 
first local maximum as the peak when describing warps. The 
sharp edge has been truncated in Figure ^ and others below 
to best show this first maximum. We will see below, that de- 
pending on the halo structure, both larger and smaller warp 
distortions are possible with the same satellite perturbation. 

4.4 In-plane response of the disk 

To complete the example, we describe the concurrent in- 
plane distortion for m = 1 and m = 2 harmonics. The overall 
distortion is small in the inner disk, Rg ^ Ro = 8.5 kpc, and 
dominated by the m = 1 term with a relative density ampli- 
tude of roughly 1.6% (4%) for the small (large) LMC mass 
estimate. The m = 1 distortion appreciable in the outer 
disk, reaching 16% (40%) near Rg > 16 kpc and increasing 
beyond that. This distortion will produce an observably lop- 
sided disk outer disk (see Figs. |ll|-[l2|). Unfortunately, it is 
difficult to determine precise distances to gas in the outer 
galaxy and this signature will be difficult to affirm. 

The quadrupole leads to a measurable oval distortion 
only for Rg > 20 kpc of 6% (16%) and is at the percent level 
or smaller near the solar circle (see Figs. [l3|-|l4[). This mild 
oval distortion is much smaller than and will be swamped 
by the predicted m = 1 signature. 



Figure 12. As in Fig. |ll| but combined with the background for 
the larger LMC mass estimate (right). 




Figure 13. As in Fig. ^ by for the m = 2 density distortion. 
5 DISCUSSION 

There are no simple formulae describing the warp ampli- 
tudes in general because of the complexity of the interac- 
tion. Qualitative guidelines are as follows. Within the range 
of scenarios explored here, the most important condition is 
the coincidence of wake pattern speed and the disk bending- 
mode frequency. Exploration suggests the 2:1 resonance be- 
tween the satellite and orbital azimuthal frequencies and the 
ILR-like resonance are most important. A secondary consid- 
eration is the location and amplitude of the wake itself which 
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Figure 14. As in Fig. |ll| by for the m = 2 density distortion. 
The mass of the LMC is exaggerated by a factor of 16 to illustrate 
the shape of the feature. 



depends on the halo profile and the existence of low-order 
resonances in the vicinity of the disk. These two features 
are not independent. However, if one could fix the orbital 
frequencies of halo stars, the wake amplitude would be pro- 
portional to the halo density and if one could fix the density, 
the wake location would scale with the orbital frequencies. 
The halos considered here are chosen to have flat rotation 
curves between the outer disk and satellite pericenter. Be- 
cause the wake is dominated by the 2:1 resonance, the wake 
peaks at roughly half the pericenter distance. Therefore, in- 
creasing the mass of the halo will tend to increase the wake 
amplitude but can decrease the disk warp if the frequency 
match with the global disk modes is less favorable. 

Because of this complicated interplay and sensitivity to 
the actual disk and halo profiles, I will illustrate the range 
of possibilities with some examples rather than give an ex- 
haustive set of models. 



5.1 Disk model dependence 

For a standard halo and satellite interaction, the disk warp 
in a Hunter A'^ = 16 and Hunter- Toomre 16X provide a 
telling comparison. The Hunter disk has a shallower less 
centrally concentrated profile than the Hunter- Toomre disk, 
naively suggesting that the Hunter disk will be more sus- 
ceptible to an outer disturbance. But, because the lower- 
frequency mode in the Hunter- Toomre disk better couples 
to forcing by the halo wake, its response is larger. Figure 
shows the response for the Hunter disk for comparison with 
Figure Figures ^ and ^ compare the same scenario for 
a mo re centrally concentrated halo model (King Wo = 7, cf. 
§5^). As above, the Hunter = 16 is less warped than the 
Hunter- Toomre 16X disk. 

Figure [is] contrasts the the Hunter- Toomre 16X and the 
Hunter A*' = 16 profiles. The latter is scaled to Rdisk ~ 3 
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Figure 15. Warp height for the Hunter TV = 16 disk with outer 
radius Rdiak = 4 and halo model Wo = 3,R = 28, M = 10 (cf. 
Fig. HI. 
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Figure 16. Warp height for the Hunter N = 16 disk with outer 
radius Rdisk = 3 and halo model Wo = 3, R = 28, M = 10 (cf. 
Fig. |l^. 



midplane height 
0.05 
0.03 
0.01 
-0.01 
-0.03 
-a05 




Figure 17. Warp height for the Hunter TV = 
radius Rdisk = 4 and halo model Wo = 7, R ■ 



16 disk with outer 
--28,M = 10. 
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Figure 18. Same as Fig. ^but for tlie Hunter- Toomre 16X disk. 




Figure 19. Comparison of the 16X disk (solid) and tiie Hunter 
A'^ = 16 disks for Rdisk =3,4 (long dash, short dash). 



and 4. Note that the Hunter A'^ = 16 for Rdisk = 3 is quite 
similar to the Hunter- Toomre 16X with Rdiak = 4 inside of 
7? ~ 2.8. Figure |l^ shows that the warp in a Hunter A'^ = 16 
disk with the same standard halo satellite perturbation is 
much smaller! Although the disk less extended, the mode 
location and morphology are to blame. The overall weaker 
halo support means that the tipping mode is closer to a bod- 
ily tip and has very low frequency. Therefore, it poorly cou- 
ples to the wake. The first retrograde and prograde modes, 
which are warp-like, have relatively high frequencies which 
also couple poorly. 



5.2 Effect of halo mass and satellite orbit 

Figures |2o| and |2l] show the warp in the Hunter- Toomre 
16X disk for increasing halo mass, M = 10, 15 and 20 re- 
spectively, for the standard model. As the halo mass in- 



creases, the amplitude of warp amplitude decreases due to a 
more poorly matched halo-wake pattern speed. In a related 
scenario, the amplitude of the wake increases with decreas- 
ing satellite energy for a fixed halo and disk as expected. 
However, the warp shape changes with energy because the 
increasing orbital frequencies couple differently to the con- 
tinuum disk response. 



5.3 Dependence on halo profile 

Here, we explore a more concentrated halo, a Wo = 7 King 
model, and therefore a sub- maximal disk model. The Wo ~ 7 
halo profile is too centrally concentrated to be an acceptable 
match to standard rotation curves: the rotation curve rises 
too steeply and over-supports the disk in the inner galaxy. 
The response for both the 16X and A*' = 16 modes are shown 
in Figures 0and|l| for both disks. The warp has a feature 
inside of of J? = 1 (inside the solar circle in the Milky Way) 
due to the larger effect of the halo at smaller radii. The outer 
warp profile is more gradual and relatively larger at smaller 
radii. Similar to the standard model, the warp in Hunter- 
Toomre 16X disk is stronger than in the Hunter A*' = 16 
disk. 



6 N-BODY SIMULATIONS 

The semi-analytic method of §^ is ideally suited to N-body 
simulation. The same biorthogonal bases can be used to rep- 
resent the potential and force fields of a particle distribution 
(e.g. Clutton-Brock 1972, 1973, Hernquist & Ostriker 1992). 
Because the Poisson equation is linear, each component can 
be represented by an expansion suited to its geometry. Sim- 
ilarly, one can easily tailor the inter-component forces to 
isolate any particular interaction. In particular, the n-body 
tests are not hampered by the difficultly in producing a per- 
fect time-dependent equilibrium disk-halo model; the disk 
feels the background halo but halo does not react to the 
background disk, as described in 



3.2 



There are two differences between the simulation and 
the perturbation calculation. First, the disk must be thick- 
ened to keep it stable against local instability (Toomre 
1964). We use a simple isothermal disk with (Tz = 20kms~^. 
This disk would also be bar unstable without a bulge com- 
ponent. A rigid Hernquist (1990) model with scale length of 
1.4 kpc and mass of 0.2 Mdisk (roug hly 1.2 X 10^° Mq) sup- 
pressed bar growth. Second, the direct acceleration by the 
satellite differentially accelerates the disk causing the expan- 
sion center to drift from the center of mass. This technical 
difficultly was circumvented for the purpose of these tests by 
dividing the original satellite in two and placing half of of 
the mass an mirror in the same orbit but 180° out of phase. 

Simulations with 100,000 particles with several different 
partitions between the disk and halo components were per- 
formed on a network parallel n-body code using the LAM 
implementation (Burns, Daoud & Vaigl 1994) of MPI (e.g. 
Gropp, Lusk & Skjelum 1994). This particular force evalua- 
tion scheme lends itself to workstation clusters because the 
communication overhead is very low (only coefficients need 
to be passed, e.g. Hernquist, Sigurdsson & Bryan 1995) as 
long as the load per node remains balanced. To suppress 
transients, the satellite is slowly turned on over an orbital 
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Figure 20. Warp height for 16X disk with outer radius Rdisk = 4 and halo model Wo = 3,R = 28, M = 15. Scaled to Milky Way, 0.1 
units corresponds to 175 pc for the large LMC mass estimate (cf. Fig. ^for M = 10). 



Figure 21. Same as Fig. Fig. ^ but for M ■ 



20. 



period. In general, the predicted features were observed: ver- 
tical wakes with midplane heights of roughly 500 pc in the 
outer disk in the predicted orientation. 

Unfortunately even with 100,000 particles and the ap- 
proximations above designed to clearly isolate the m = 1 
wake, the simulation results are difHcult to interpret due 
to discreteness noise. Simulations without any satellite re- 
veal that the noise has two effects. First, the disk origin 
random walks in the fluctuating Z = m = halo compo- 
nent. This causes m = vertical distortions which lead to 
scale height thickening. One should expect these low-order 
fluctuations are amplified above the Poisson amplitude by 
the global gravitational response as described in Weinberg 
(1993, 1997). Secondly, one finds m = 1 height distortion 
from the I = 2 noise-excited halo component similar in scale 
to that expected from the satellite excitation. In retrospect, 
this might have been expected. Any response to a pertur- 
bation will feature any discrete modes. For the halo, the 
strongest of these will be weakly damped. The halo wake, 
which is now the self-gravitating response to the disk distor- 
tion, will be contain the same weakly damped mode found 
in the response to the satellite. Therefore, particle fluctua- 
tion noise can excite a similar with modes and a similar disk 
response. The amplitude of the noise excited distortions was 
as much as 30% of the amplitude observed with a satellite 
present. However, different Monte Carlo realizations of the 
same initial model gave vertical warps of arbitrary orienta- 
tion whereas in the presence of the satellite, the predicted 
orientation is obtained. 

A proper n-body investigation of these effects will re- 
quire simulations with many more particles and will be left 
for the future. This investigation, however, inadvertently 
raises the interesting possibility that a satellite may interact 
with noise-excited features to produce an overall response 
which is larger than the time-asymptotic predictions de- 
scribed above. Because the modes are already excited by 
noise, they may be entrained and strengthed by the exter- 
nal perturbation. This is similar to other noise amplification 
mechanisms found in Nature, e.g. stochastic resonance mech- 
anism observed in neurobiology (e.g. Collins, Chow & Imhoff 



1995). Moreover, there are many sources of inhomogeneity 
in real galaxies such as star clusters, gas clouds, unmixed 
streams from dwarf dissolution, and other sources of per- 
turbations such as tidal distortions from interactions with a 
host galaxy cluster and near-neighbor interactions. In fact, 
the noise in the simulation presented here corresponds to 
that produced by a halo of IO^Mq black holes (e.g. Lacey 
& Ostriker 1985). Altogether, it seems likely that observable 
effects on disks due to halo interactions will be stronger than 
the predictions in §kl 



7 SUMMARY AND FUTURE WORK 

This paper summarizes the effect of a satellite companion on 
a galaxian disk embedded in a responsive (or live) halo. Be- 
cause the mechanisms described here involve multiple time 
and length scales, the interaction is complex and depends on 
the details of each dynamical system: satellite orbit, halo and 
disk structure. The general conclusions and expectations are 
enumerated below: 

(i) The halo can sustain a significant wake in the presence 
of an LMC-like satellite. This wake is most often caused by 
a 2:1 orbital resonance and therefore peaks roughly half way 
between the satellite orbit in the halo center. 

(ii) The halo wake, because it has structure at smaller 
galactocentric radii, can excite warp modes in the disk more 
efficiently than direct tidal perturbation. 

(iii) The strongest warps obtain for a near commensura- 
bility between a disk mode and the pattern speed of the halo 
wake. This makes robust predictions for warp amplitudes 
more difficult, but we find that warps with observable am- 
plitude can easily result by the mechanism described here. 

(iv) For halos with roughly flat rotation curves out to 
roughly 50 kpc, the satellite orbit (pericenter, plane orienta- 
tion, and eccentricity) which determines the forcing frequen- 
cies and disk profile which determines the bending frequen- 
cies are most important in determining the induced warp. 

(v) A polar satellite orbit will produce the largest warp. 
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The inferred LMC orbit is nearly optimal for maximum warp 
production. 

The halo wake (e.g. Figs. |and||) plays a critical role in pro- 
ducing a warp and is an observable consequence of the mas- 
sive dark halo. A warp survey combined with recent satellite 
surveys (Zaritsky et al. 1993) may provide additional statis- 
tical evidence for the massive halo hypothesis. 

Several possibly important interactions have been ig- 
nored in this study and present topics for further research. 
First, the disk is model is infinitely thin; including the three- 
dimensional structure will damp the warp producing modes 
but this damping is likely to be small at large scales (Wein- 
berg 1991). Second, the disk respond to the halo wake and 
the halo responds to the two-dimensional disk distortion 
but not the three-dimensional one. Therefore the dynami- 
cal friction against the halo explored by Nelson & Tremaine 
(1995) is not included. The calculational method used here is 
tractable because of the assumption that all transients have 
mixed away. As they mix, they produce fiuctuations on many 
scales. Intrinsic halo inhomogeneities such as clouds, clusters 
and massive black holes are also a source of noise. Together, 
this noise may seed interesting observable features, such as 
inner bars and arms, which are not part of the long-term 
wake. In addition, as described in these fluctuations will 
drive the same modes which produce the wakes at largest 
scales. N-body simulations suggest that noise-excited struc- 
ture can have a effect on the disk similar to the satellite 
companion. This leads to the possibility that the amplitude 
of the large-scale response to an external disturbance may 
be amplified by entraining the pre-existing noise-excited fea- 
tures. 
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APPENDIX A: COMPUTING THE DISK 
DISTRIBUTION FUNCTION 

The distribution function is the solution to the following 
integral equation: 
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= / d\f{E, J), 



(Al) 



where E and J are the orbital energy and angular momen- 
tum, respectively. For a given energy, the maximum angular 
momentum, that of a circular orbit with energy E, is de- 
noted by Jmax{E). The potential of the halo {^haio) and disk 
profile (E(_R), ^^disk) are fixed and assumed to be known to 
start. We allow the distribution function to be represented 
by a Gaussian basis in the variables E and k= J / Jmax (E) : 

f{E,K) = Y,a,,e^p[-{E-Eif/2al-{K-ii,f/2al] 



(Bl) 



Because r(f) is quasiperiodic in two frequencies for a spheri- 
cal halo, we may expand equation (Bl) (or eq. ^in the main 
text) in a Fourier series in time. 

This is straightforwardly done in the perturber's orbital 
plane and then rotated to the desired orientation using the 
rotational properties of spherical harmonics (e.g. Edmonds 
1960): 



(B2) 



(A2) 



At any point Rh, the distribution function is related to 
E(i?fc) through equation (A_l) 



E(i?fe 



dVr 



dvtfiE,K) (A3) 



where Vmax{Rk) = ^/^{Emax 



-I'd 



^halo), E = {vi.+ 



Vt)/2 + <bdisk + ^haio, and k = RkVt/Jmax{E). An exact 
solution requires that the aij to satisfy equation (Al) for 



all R. We can state this demand as the set Uij which min- 
imizes the square of the difference of T.{Rk) and E(-Rfc) at 
all R. This demand may be discretized to the set Uij which 
minimizes 



X 



= ^Wfc [E(i?fc)-S(7?ft)]' 



E 



■Wk 



'^^aijT,ij{Rk) — t,{Rk 



(A4) 



with the condition 

Ke,k)>0 



(A5) 

where Wk is a weighting, which may be Wk ~ 1. This de- 
fines a standard quadratic programming problem for the aij 
which we solve using Powell's algorithm (1982) as imple- 
mented in the QLD code provided by Andre Tits (Lawrence, 
Zhou & Tits 1994). 

Distribution functions used for models described here 
used a grid of 20 x 2 in E and n and penalized the expres- 
sion in equation ( A4 ) to construct a tangentially anisotropic 
distribution, 



2 



E 



■Wk 



E' 



'ij ^ij 



i,j r,s ) 



+ 



(A6) 



with itJfc = 1, A = 10 and a — &. 



APPENDIX B: PERTURBATION 
COEFFICIENTS FOR AN ORBITING 
SATELLITE 

The biorthogonal expansion coefficients (cf. eq. |^ for a satel- 
lite orbiting in a spherical halo are conveniently derived by 
expanding its gravitational potential in an action-angle se- 
ries. Assuming a point-mass perturber, the coefficients are 



where 



(B3) 



The r\^f.{(3) are the rotation matrices and a, f3, 7 are the 
Euler angles describing the orientation to the orbital plane. 

Using this and the inverse of equation (B2), we can now 
expand equation ( Bl ) in action-angle variables. For spherical 
systems, the third angle variable describes the line of nodes; 
it has zero frequency and has been suppressed. The action- 
angle coefficients are then defined by 



Oil 



(2^)2 ^ 



e'^"Yik{n/2,0)e'^^pr'ir{w^)) 

— J2Yk{7v/2,0)rU{l3)5ki, {l3)e X 

k 



dwie ^ Pi [f'KWi))- 



(B4) 



The angles wi and W2 describe the equal-time motion from 
pericenter to pericenter and the mean azimuthal motion re- 
spectively. The angle a describes the rotation of the orbital 
plane; q = corresponds to the y axis in the orbital plane 
coincident with the line of nodes. The angle 7 describes the 
orientation of the line of nodes in the original azimuthal co- 
ordinate; 7 = places the line of nodes coincident with the 
y axis in the original system. Because the orbital radius and 
the deviation of the true azimuth from the mean azimuth 
(-0 — ii)2) is even in wi, the integral over w-i may be simplified 
and one finds: 



67 = e''"^e"^°y„,(^/2,0)rU(/3)H//?;;(I), 
where W is defined by 



(B5) 



Wli'il) = - dwi cosihwi + l2f{wi))p'rir{wi)) (B6) 
^ Jo 

with 

f{wi) ^ jdr [2{E - $0) - j'/r')] {0,2 - J/r") . (B7) 
Each term in the expansion in equation (|^) is then 



irwy il2(y. 



ll ,l2 — — oci 



e"'"Yi,4Tr/2,0)W;i 



(B8) 



for wio = W20 = 0. We will focus on the m = 1, 2 responses 
for I < 4. Several tests with I < 6 suggest that I < 4 domi- 
nate the large-scale response. 
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APPENDIX C: NUMERICAL COMPUTATION 
OF RESPONSE MATRICES 



The elements of the response matrix have been derived in 
several places (Paper II and references) and are given by: 

Tjim, wo \3 2 Y^Y^ [ dEdJJ dfo 

7^., (.) = ^^1^1^ J TME;Tf'-^ ^ 

j^^\Yn. (V2, 0) r (IX,^ (I). (CI) 

The inverse Laplace transform as described in §^ requires 
evaluation of matrix elements of this form in two contexts. 
The Laplace transform of the coefflcients in equation (|^, 
6'™'(s), contribute a simple pole on the real axis is easily 
performed by deforming the contour to 9(s) ^ —<x. The 
inverse transform then has the following form 

^ ' 2-Ki 



Tvi I 

= 77- / dse"*0~'(s)7^(s)■b''" — 
2tvi / - s — 

J c — zoo 



.(C2) 



Although the stabihty assumption ensures that has no 
poles in the 3ff(s) > half plane, there are poles in the 
5R(s) < half plane corresponding to damped modes (e.g. 
Weinberg 1994). Of course, weakly damped modes will be 
subdominant to an oscillatory mode after sufficiently long 
time and we assume this limit here. 

The elements of T>{s) for the deformed integration path 
must be analytically continued to 5R(s) < from 5R(s) > 
where the transform is defined. This leads to a Cauchy inte- 
gral with two simple poles on the real axis for each matrix 
element (cf. eq. |Cl[) : s — icj and s = il ■ ft. Numerically, 
this may be straightforwardly evaluated by subtracting the 
singularity from the integrand and evaluating it separately. 
Consider the following Cauchy integral: 

dzJ^ = / dzM^I^ + f[zo) I dz^. (C3) 



The first term on the right hand side is non-singular and can 
be evaluated by simple quadrature. The second term can 
be performed analytically once the contour is specified. In 
this case, analytic continuation requires the standard Lan- 
dau contour (cf. Krall & Trivelpiece 1973). Unlike the case 
of damped modes, the evaluation of the forced response re- 
quires no extrapolation or explicit complex integration. 

This paper has been produced using the Royal Astronomical 
Society /Blackwell Science style file. 
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